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Compressible Seal Flow Analysis Using the Finite 
Element Method With Galerkin Solution Technique 

by 

John Zuk 

NASA -Lewis Research Center 

ABSTRACT 

High pressure gas sealing Involves not only balancing the viscous 
force with the pressure gradient force but also accounting for fluid 
inertia — especially for choked flow. The conventional finite element 
method which uses a Raylelgh-Ritz solution technique Is not convenient 
for nonlinear problems. For these problems, a finite element method 
with a Galerkin solution technique (FEMGST) was formulated. One 
example, a three-dimensional axlsymmetrlc flow formulation has non- 
linearities due to compressibility, area expansion, and convective 
inertia. Solutions agree with classical results in the limiting 
cases. The development of the choked flow velocity profile is shown. 
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Nomenclature 

2 2 

cross-sectional area, in. ; m 

basis function geometric coefficient evaluated in 
Table 2 

coefficient expression used in GADE, defined in Table 2 

basis function geometric coefficient evaluated in 
Table 1 

coefficient expression used in GADE, defined in Table 2 
undetermined coefficient 

basis function geometric coefficient evaluated in 
Table 1 


am 


S M 

"am 


pXV 

m 


unknown coefficient matrix 

coefficient expression used in GADE, defined in Table 2 

coefficient expression used in GADE, defined in Table 2 

linear basis function 

basis function of an element m 


XY 

F coefficient expression used in GADE, defined in Table 2 

-am * 

XY 

G coefficient expression used in GADE, defined in Table 2 

am 

h film thickness (gap), in. ; m 

N total number of mesh points 

n an integer 

2 

P static pressure, psi; N/m 

R basis function integral defined in Table 4 

basis function integral defined in Table 4 

R^ basis function integral defined in Table 4 

RZ basis function integral defined in Table 4 

RZ^ basis function integral defined in Table 4 


Jt 



leakage flow Reynolds number in radial direction, pUh/P 


modified Reynolds number, Re(h/AR) 
radial direction coordinate 
dimensionless radial coordinate, r/AR 

temperature, °F; K 
dimensionless temperature, T/AT ^ 

leakage flow reference velocity, ft /sec; n/sec 

velocity in r-direction or x* direction, ft/sec; m/sec 

dimensionless velocity, u/U 

dependent variable identification number 

number of dependent variable 

coordinate in pressure gradient direction 
basis function integral defined in Table 4 

basis function integral defined in Table 4 

basis function integral defined in Table ■+ 

basis function integral defined in Table 4 

coordinate across film thickness 
dimensionless coordinate, z/h 

dependent variable 
basis function 

area of triangle 

error, global or series trunkation 

absolute or dynamic viscosity, (lbf) (sec) /in.^ ; 

(N) (sec) /nr 

2 2 

kinematic viscosity, ft /sec; m /sec 
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P density, (lbf ) (sec 2 ) /in. 4 ; kg/m 3 

P* dimensionless density, p /p Q 

2 o 

T shear stress, lbf/in. ; N/cm 

{ i Fn ^ basis function integral defined in Table 3 
Subscripts: 

i vertex number of triangle element 

j vertex number of triangle element or dependent variable number 

k vertex number of triangle element 

L triangle element identification 

m mean or triangle sub-element region 

N location along flow leakage length 

n vertex number 

U triangle element identification 

v variable number 

x vertex location coordinate, axial direction 

y vertex location coordinate, radial direction 

a vertex identification of triangle element 

6 vertex identification of triangle element 

Superscripts : 

X vertex location coordinate, axial direction 

Y vertex location coordinate, radial direction 

* dimensionless quantity 
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INTRODUCTION 

A purpose of chis paper Is to devel p a numerical analysis tool 
which will enable analysis of complex seal problems. The numerical 
method of solution utilizes a new procedure — the finite element method 
with a Galerkin solution technique (FEMGST) rather than the method of 
finite differences. Utilization of this FEMGST should enable the fluid 
film analysis to be readily integrated with solid mechanics and dynamics 
finite element analyses which are carried out by existing computer pro- 
grams. Also, the Galerkin solution procedure can solve nonlinear 
problems which cannot be solved using the conventional finite element 
method which uses a Rayleigh-Ritz solution procedure (1) . 

In pressure balanced face seals and self-acting lift pad seals (2) 
(Figs. 1 and 2), the flow can choke. Thus fluid Inertia must be ac- 
counted for In a flow analysis. Since inclusion of fluid inertia makes 
the flow equations nonlinear, the flow is currently solved using an 
approximate integrated average model (3). 

These approximate seal analyses may not yield satisfactory results 
when the nonlinearity becomes large. Further, these approximate methods 
also do not yield detailed distribution information about the flow field. 
Thus it is necessary to utilize a numerical solution scheme. A finite 
element method with a Galerkin solution technique is used. A finite 
difference scheme could have been used. However, it was felt that 
ultimately a finite element seal flow analysis would be a powerful 
analysis tool. 
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Some of the advantage s of the finite element method over the finite 
difference method are: 

1. Boundary conditions need not coincide with coordinate lines. 

2. Layout of nodal points is more flexible. 

3. Boundary conditions of abrupt geometric changes, such as step 
changes in film thickness, do not require finding a special inter facial 
boundary condition such as mass flux continuity (4). 

4. The method can yield numerical solutions more rapid) y. Studies 
on a transient one-dimensional heat conduction problem show that this 
method is about six to 60 times faster than the finite difference method, 
depending on the desired accuracy (5). 

The disadvantages of using this method are: 

1. The solution is an average value over the element rather than 
an exact solution at every point in the region; however, this is 
satisfactory for most engineering problems. Better accuracy can be 
obtained by decreasing the element size. 

2. Programming is also more difficult. A larger software effort 
may be required. 

3. The coefficient matrix for the nonlinear flow problem is non- 
diagonally dominant due to the momentum equation. 

Although seal flow problems, in particular, have not been analyzed 
using finite element techniques, other lubrication flows have been 
analyzed using the finite element method. Reddi (4) and (6) has solved 
some incompressible and compressible bearing lubrication problems using 
the finite element method with the Rayleigh-Ritz approximate solution 
procedure. The Reynolds lubrication equation in quadratic functional 


y 
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form and thus is an admissible form to a Rayleigh-Ritz solution. 

Recently Oh and Heubner (7) have applied the finite-element techniques 
to solve the elastohydrodynamic finite journal bearing problem. The 
Reynolds' equation for the fluid film and the three-dimensional elasticity 
equations for the bearing housing were solved simultaneously using an 
Iteration scheme. 

The usual finite element method, which is based on a Rayleigh-Ritz 
approximation, has been extended to other continuum problems such as 
heat conduction and dynamics (1). This means that a variational principle 
must be found; e.g., the total potential energy of a system is stationary. 
However, in fluid mechanics such a principle may not be readily found 
in a workable form. Also, the classes of nonlinear problems that can be 
solved are very limited. The problem must be reducible to a quadratic 
functional. Such problems are described in equations which are sometimes 
called "quasi-harmonic" equations. The heat conduction equation with 
variable conductivity, irrotational flow of ideal fluids, Lagrange's 
equation, and the compressible Reynolds lubrication equation are examples 
of "quasi-harmonic" equations. The compressible Reynolds' equation, 
however, neglects convective inertia forces whose retention along with 
the viscous forces changes the character of the equation so that it is 
no longer "quasi-harmonic". 

The Galerkin method (8) by itself could be used. Cheng and Pan (9) 
applied Galerkin' s method in solving the nonlinear unsteady Reynolds 
equation. It allowed the reduction of Reynolds equation directly from 
a partial differencial equation to a system of first-order, ordinary, 
differential equations, which together with the equations of motion 
of the journal bearing, yielded a tractable stability analysis of 
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finite plain journal bearings. Unfortunately, in many problems the 
series is unworkable when the variables are changing rapidly. Many 
terms are required in order to have a satisfactory solution. This 
problem may be overcome by using the Galerkin method with piecewise 
continuous derivative basis functions. Thus the formulation in this 
paper is essentially the Galerkin method with finite elements as the 


basis functions. 
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RESULTS AMD DISCUSSION 

The finite element method with Galerkln solution procedure will be used to 
solve 3-dlmensional axisymmetric and 2-dimensional seal flow problems. Ar. 
overview of the technique will be presented by outlining the general procedure. 
The method is basically as follows: 

1. Subdivide the region into triangular ring elements as shown 

in Fig. 3 . 

2. Define an "irregular pyramid" basis function F n (x,y) at 

each nodal point. (A linear basis function - first order splines.) 

See Fig. 4. Within each triangular element, F n is a linear 

function of distance such that F =1 at n * (x,y) and F * 0 

n n 

on the boundary. Outside of the region defined by the triangular 

elements touching the point (x,y), * 0 . 

Where: n is the nodal point number, n = l,***, N; (x,y) 

describe the nodal points in terms of their radial 

and axial locations. That is, x is the radial 

direction nodal point number, x = 1,***, X * M, 

while y is the axial direction nodal point num- 

ber, y = 1,*“, Y - N 
max 

3. Let 

N 

\ \ 

fj(x,y) = ^ / C nj F n (x,y) M 

n*l 

Where: j is the dependent variable number, j - l,***, J 

4. Place the piecewise continuous series into the governing 
equations. 

L^jCx.y)} = [ 2 ] 

Where: i is the equation number, i * l,***, J 

5. Orthogonalize the global error to each basis function as re- 


quired by the Galerkin method. 



10 


6. Place the boundary conditions on the appropriate nodal 
points . 

7. Solve the resulting NJ linear cr nonlinear algcbr-a' eq* 
tions for NJ unknowns. 


The details of the finite element method with Galerkin solution 
technique will now be discussed. The domain of interest is subdivi- 
ded into triangular subregions as seen in Fig. 3 (a plan view of 
the r-z plane). Consider a triangular subregion with vertices i, 
j, and k as seen in Fig. 5 * 

Let the dependent variables » ^ 2 * ^3* * * * ^^ near ^V 

varying in this triangular subregion; i.e., let 


* 

* 


1 

2 


a l + ^l r + C 1 Z 


a 2 + V + C 2 Z 


a V + b V r + C V Z 


[3] 


The linearly varying field assures continuity between elements 
since ’’lines which are initially straight remain straight in their 
displaced position”. The nodal or vertex values of ^ are defined 
uniquely and continuously throughout the region; however, the first 
derivative is discontinuous, hence, this function with piecewise 
continuous derivatives is a first order spline. At the nodal points 
or any triangular element, a dependent variable can be expressed as 
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* 

* 

n 


li 

U 

lk 


»1 + b l r i + C l Z i 


a l + Vj + Vj 


a l + b l r k + C l Z k 


[A] 


which in matrix form becomes 

{*. } » [D] {a} 15] 

L n 

or 

{a} - [D] I {» 1 } [6] 

1 n 

whe re 



Hence, the unknown coefficients from Eq. [1] have been expressed 
in terms of the vertex position and dependent variable values at 
these nodal (vertex) points. The inversion of the coefficient ma- 
trix, [D], is shown In (10). Hence, it can be readily shown 
that the resulting form of the dependent variable in a triangular 
subregion m is: 


v (m) 


1 

2A 


[o 


V 


C i z) *vi 


+ (a, 


b.r 

1 


c j 2) Vj 


+ (a k + V + c k z) *vk] 181 

m 


v 

i 



The values for a, b, and c are shown in Table !• The basis 
function can also be expressed in the above form. The basis func- 
tion equation is developed in (10) and is 
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F u " ie (a i + b i r + c i z)F i + (a j + y + c j z)F j 


+ ( \ + b k r + c k z)F k 


[9] 


where the basis function may have one of the following combinations: 

{F i -1, Fj * 0 , F k » 0), {F i »0, Fj * 1 , F k * 0}, * 0, F - 0, F fc - 1} [10] 


A given basis function whose value is unity at nodal point XY, will 
have an irregular pyramid’s base as its region of definition. The 
convention that is used results in six member triangular subele- 
ments - which comprise the region of definition of a basis function 
located at an interior mesh point. See Fig. 5. 
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As previously described, the piecewise continuous derivative 
representation of each dependent variable, 18], is placed Into the 
governing Eq. [2] 


t 


i 


X 
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Then the global error is orthogonalized to each basis function as 
required by the Galerkin method, Eq. [3] 


/ 


e,F dA * 0 
i n 


Since the basis functions, F , are zero outside of the irregu- 

n 

lar base pyramid region, Eq. [3] can be written as 


M,N 

£ 

n=(X,Y)=(l,l) 


e.F dA = 0 
r n 


[123 


Thus, a term of the sum in [12] becomes the following for the 
basis function that has a value of unity at n * X,Y, the apex of 
the pyramid. 


e.F dA + /e.F dA + / e.FdA+ / e.F dA 

i xy / x x y / i xy / i xy 


' L-l 


L+l 


U-l 


+ /e.F dA + / e.F dA * 0 

l xy / l xy 


[13] 


'U 


'U+l 


At a point of any triangular region of the pyramid base, a 
dependent variable is influenced by the values at its vertices. 
This can be expressed as 


* = c i F i + Vj + c k F k 


114] 


V* • * 
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Since etc., Eq. [14] can be expressed as 


* " Vi + Yj + Vk 


where: F , F , F, are basis functions that have the apex * unity 
i j k 

at vertices i, j, k, respectively; 4^, 4^, 4> k are the magnitudes 
of the basis vectors which describe the plane (equal to the magni- 
tude of the variable at the vertex) . 


Sy utilizing Hue a a. basis fuuk.uiuu&, cue icaullaul ^oluciim 


vector (where the coefficients have been determined) describing the 
triangular plane will remain a plane; at the interface of the tri- 
angular element (boundaries) the functions are continuous but not 
the derivatives (since the basis function is a first order spline). 
The whole volume of each basis function can be geometrically inter- 
preted as an irregular hexagonal based pyramid function. 

See Fig. 5 for a plan view of the region of definition of a typ- 
ical basis function, f^. Fig. 6 illustrates a 3-dimensional geo- 
metric interpretation of the function space. Note, the bases are 
irregular to show a variable mesh situation. 



SEAL FLOW EQUATION FORMULATION 


The compressible seal leakage flow will not be analyzed using the 
finite element method with Galerkin solution technique (FEMGST). 

The proper dimensionless governing equations where "quasi-fully 
developed'*, parallel flow exists, using shear stress as a dependent 
variable, are: 

1. Conservation of mass 

jf* <p*r*u*) - 0 [15] 


2. Conservation of momentum 


Re * p * u * is* _ bt* . o 

Ke p u 9r* 3r* 3z* u 


3. Perfect gas relation 


p *t* - p* * 0 


4. Shear stress-velocity relation 


3u* 

3z* 


— T* * 0 


where 


Re* * Re 


(v R l) 


Only Isothermal flow cases will be considered. 


[16] 


[17] 


[ 18 ] 


FEMGST Form of Equations 

The procedure for forming the Galerkin approximation equations for 
the above set of partial differential equations will not be presented 
using the procedure described in the previous section. Details of this 
procedure can be found in (10). 
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The equations governing isothermal compressible leakage flow in 
FEMGST form at nodal point XY are 
1. Conservation of mass 


m , 


am 


m=L-l 


b y 

2b 

m 


p u 
Y a 


+ < - + 

a3m a 3 


3 xy 

am 



= 0 


[19] 


2. Conservation of momentum 


U+l 


r « * (\ \ . XY„ _XY 

Re Itt-I u B „ u p 0 + E P_ - F 


J 


m=L-l 


2b y a3m a 3 am a ^*m a 
, m ' 


- 0 


3- Shear stress-velocity relation 
U+l 


L / \ram a am a/ 


m=L-l 


[ 20 ] 


[ 21 ] 


4. Equation of state 
U+l 

T' (b*? - 0 [22] 

J j \ aSm a 3 am a/ 1 J 

m=L-l 

_ , a XY ^XY T{ w XY ^XY j „XY . , 

The functions A , B 0 , D , E , F , and G are defined 
am a3m am am am am 


in Tables 2, 3 and 4. 

A typical seal flow region subdivided into triangular ring elements 
is shown in Fig. 3. The flow is treated as axisymmetrlc. Basically, 
FEMGST is the Galerkin method with finite elements as basis functions. 

As previously mentioned, at each nodal point an "Irregular pyramid" 
basis function F (r,z) is defined (see Fig. 4). 


v 

i 


v- 
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DISCUSSION 

As a starting point, the first case studied in detail was incom- 
pressible radial Poiseuille flow. This flow is simpler in nature than 
the cases of interest. For this case Eqs, [15], [16], and [18] are 
solved with Re* - 0 and P* - 1. 

Figure 7 shows the model and conditions used for solving incom- 
pressible radial Poiseuille flow. Figures 8, 9* and 10 compare the 
numerical solutions with the exact analytical solutions at various 
radial locations for velocity, pressure, and shear stress, respectively. 
Note, there is excellent agreement for the mesh size chosen. 

In order to check the compressible flow formulation of FEMGST, 
a special inviscid flow problem was solved. The modified Reynolds 
number Re* was a parameter varied in this study. The inviscid flow 
problem had all of the salient features of the general compressible 
flow seal leakage problem. That is, nonlinearity due to compressibility 
area expansion, and convective inertia. However, the problem is simpli- 
fied because the shear stress is everywhere zero. The radial velocity 
distribution is shown in Fig. 11. There is excellent agreement between 
the exact and FEMGST solutions. 

A compressible seal flow case of practical interest was then solved 
A seal was studied with a radius ratio of 0.8 (see Fig. 12). The 
transverse velocity profiles at the seal exit are shown in Fig. 13, 
whereas the radial pressure distributions at the dimensionless trans- 
verse gap value of 0.25 are shown in Fig. 14. The modified Reynolds 
number again was a parameter, but its value here is a measure of 
the convective inertia effect in the flow. For Re* * 0 the results can 
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be compared with the exact viscous compressible flow solution* Note 
the good agreement of FEMGST with the exact solution. For Re* + 0 there 
is no exact solution for comparison; however, the results are physically 
expected. The solutions were terminated at Re* * 12. At Re* * 12, the 
results indicate that the limiting exit sonic velocity conditions had 
been reached at the centerline. The solution had to be terminated 
because results indicated that the formulated set of equations was no 
longer valid. The radial pressure distribution in Fig. 14 shows that 
the exit pressures change the most with Reynolds number. This is expected; 
since the velocity increase is the greatest in the exit region, the 
density (pressure for this isothermal case) must correspondingly de- 
crease. This can be seen by examining the compressibility terms in 
the mass conservation equation. 

To use FEMGST, a new computer technique had to be developed. The 
system of equations from FEMGST could not be solved by the usual 
iterative solution methods. Thus, a direct method algorithm was 
developed and used to reduce the large, sparse, nonsymmetric band 
matrix. This method greatly reduced machine storage requirements and is 
described in (10). 


i 
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Table I. - Values of Geometric Basis Function 
Coefficients 


r j z k ~ Vj a j * Vi - Vk \ 


r i z j " r j z i 


b i a 2 j ' 2 k 


Z 1 - 2 . 
k 1 


\' z i- 


c l ■ r k. - r J 


C J * r i " r k 


c 3c r — r 
k r j r i 



Table II 


Coefficient Expressions Used in Galerkin Approximation 
Difference Form (GADE) 


Variable Form GADE Form 



Note: 
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Table III. - Basis Function Integrals For a Basis Function Whose 
Apex (Unit Value) is at Nodal Point a 


1. L-l Element; a 


j, therefore, F^ 


1 









[“jVi + v + c j z] 

[ajR + bjR 2 + c^RZ] 
[a^Z + bjRZ + c^z 2 ] 
UjRZ + bjR 2 Z + c^RZ 2 ] 

[a^R 2 + b^R 3 + c^r 2 z] 
[a^Z 2 + bjRZ 2 + c^z 3 ] 


Repeat above calculations for: 

2. L Element; except a » k # * 1 

3. L + 1 Element; except a = k, * 1 

4. U - 1 Element; except a * j , = 1 

5. U Element; except a ■ i, ■ 1 

6. l ? + 1 Element; except a * i # F i = 1 





1 


Table IV. - Values of Basis Function Integrals 


R * / / rdrdz ■ d A 

J kJ r 


* / i zdrdz * d A 


* // r2drdz ■ & + + $ * < 
- fj zh'iz + ■? * zl) * u * a 


K ’ TT r * drdZ ' S <r l*i + Vj + Vk> + “V. 


R 3 - yy r 3 drdz 


R 2 Z - y^^r 2 zdrdz 


RZ 2 - rz ^ drd2 


Z 3 - Z 1 / 1 z 3 d rd z 
1/ V 


1 f i Z i 


4 ■ If drdz ■ T 1 'j *, 


I 1 'k *k 
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z dzdr 


i[<i< 


r k ‘ r i ) + A jk (r j - f k ) + A ij r i “ f j ) 




223 3 2233 223 3 

+ 2A ki B ki (r k - r i> + Vi:/VV +w ijV r i • V 


344 34^' 344 

+ \i B ki (r k ' r i> + A ikV r i - 1 : A ii B ii , - r i- r i ) 


+ i B ki (r k • r i } + i B 5k (r j " r k ) + 


1 i 5 5l 

5 -iJ U i-Vj 


2 , . 1 1 a 3 . 2 213-2 2. 1 3 . 2 2, 

z rdzdr - - ? (r k - rj + y A jk < rj - r fc ) + j A lj (r i " r j> 




. 3 . _2 . 4 4.^ 3 . „2 , 4 4. 

+ 4 \i B kl r k ’ r i + 4 A jk B jk (r j " 


, 3 . ,2,4 4. . 1 d 3 , 5 5. . 1 .3 , 5 

4 A i j B lj r i r j ) + 5 B ki r k r i ) + 5 B jk (r j 


. 1 .3 , 5 5, _2 

+ 5 B ij (r i ' r j> “ Z R 


Table IV. - Continued 
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( |.r 2 d«dr ■ i [i A^fr’ - r B ) + ± A^r* - r B > * i *| j( r 2 - r 2 ) 


♦ i 


2 \i B ki (r k r i J + 2 A jk B jk (r j “ r k* 


+ 2 A ij B ij (r i ' r j ) + 5 B ki (r k " r i } + \ B jk (r j " r k } 


. 1 .2 , 5 5 . 

+ 5 B ij (r i " r j> 


ZR 


n 3 
r d*dr 


T \i (r £ - r t> + -f < r j - ^ + ^ << - 


r J> 


+ T B kl (r k - '!> * T K.i'l - rb 


T s ij (r i- r j>] 


5 “jk'‘j 
3 


where: 


h 1 h 1 S 


r i + r i + \ 


i . Il z j ~ V* . 

ij r i ~ r j - c k ’ 


Jk 


r 1 z k " Vi _ a i c l c k 


r j ‘ r k 


-c , 


ij 


. b x c i c j . 

r i " r j -C k 


, - . ii!si 

jk r j - r k - c i 


Vi - r i z k a i c i c k 


r k ~ r i 


-c 


z k - l i _ b i c t c k 


ki r fc - t, 


Table IV. - Continued 






FEMGST 


z 



Figure 3. - The finite element idealization of flow through a 
gas film seal. 



(a) THREE-DIMENSIONAL REPRESENTATION OF THE BASIS FUNCTION. 


(X-1MY+1) 


X(Y+1) 



(X+1IY 


X(Y-l) (X+1XY-1) 

(b) BASE OF IRREGULAR PYRAMID USED TO DEFINE BASIS 
FUNCTION. 


Figure 4. - Geometric interpretation of an 
irregular hexagonal based pyramid function. 
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Figure 6. - Geometric interpretation of function space (irregular hexagonal based pyramid 
basis functions). 
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Figure 7. - Incompressible radial Poiseuille flow sample problem. 
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Figure 8. - Comparison of FEMGST solution with exact solution for radial velocity distribution. 
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Figure 10. - Comparison of FEMGST solution wi* 1 exact 
solution for pressure distribution 
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Figure 11. - Comparison of exact inviscid flow solution with 
FEMGST for velocity distribution and Re* • 1, 1. 5, 5, and 
100 . 



Figure 12. - Domain and boundary conditions for compressible seal leakage 
problem. 
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Figure 13. - FEMGST solution for exit velocity profile for a range of Reynolds num- 
bers (0-12); 11 x 17 uniform mesh. 



Figure 14. - FEMGST solution for radial pressure dis- 
tribution at an axial location of z* ■ 0. 2500 for a nasa-l.wi, 

range of Reynolds numbers (10-12). 


